home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_8 / nelder.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.6 KB  |  121 lines  |  [MATF/MATL]
  1. function [V0,y0,dV,dy,P,Q] = nelder(Fn,V,min1,max1,epsilon,show)
  2. % [V0,y0,dV,dy,P,Q] = nelder(Fn,V,min1,max1,epsilon,show)
  3. % Nelder-Mead search for a minimum.
  4. % Fn is the function, input.
  5. % V  is the starting simplex, input.
  6. % min1 is the minimum number of iterations, input.
  7. % max1 is the maximum number of iterations, input
  8. % epsilon  is the tolerance, input;
  9. % show  if show==1 the iterations are displayed, input.
  10. % V0 is the vertex for the minimum, output.
  11. % y0 is the function value  Fn(V0), output.
  12. % dV is the size of the final simplex, output.
  13. % dy is the error bound for the minimum, output.
  14. % P  is a matrix containing the vertex iterations, output.
  15. % Q  is an array containing iterations for  F(P) , output.
  16. if nargin==5, show = 0; end
  17. [mm n] = size(V);
  18. for j =1:n+1,
  19.   Z = V(j,1:n);
  20.   Y(j) = feval(Fn,Z);
  21. end
  22. [mm lo] = min(Y);            % Order the vertices:
  23. [mm hi] = max(Y);
  24. li = hi;
  25. ho = lo;
  26. for j = 1:n+1,
  27.   if (j~=lo & j~=hi & Y(j)<=Y(li)), li=j; end
  28.   if (j~=hi & j~=lo & Y(j)>=Y(ho)), ho=j; end
  29. end                          % End of Order.
  30. cnt = 0;
  31. while (Y(hi)>Y(lo)+epsilon & cnt<max1) | cnt<min1
  32. % The main while loop has started.
  33.   S = zeros(1,1:n);          % Form the new points:
  34.   for j = 1:n+1,
  35.     S = S + V(j,1:n);
  36.   end
  37.   M = (S - V(hi,1:n))/n;     % Construct vertex M:
  38.   R = 2*M - V(hi,1:n);       % Construct vertex R:
  39.   yR = feval(Fn,R);
  40.   if (yR<Y(ho)),
  41.     if (Y(li)<yR),
  42.       V(hi,1:n) = R;         % Replace a vertex:
  43.       Y(hi) = yR;
  44.     else
  45.       E = 2*R - M;           % Construct vertex E:
  46.       yE = feval(Fn,E);
  47.       if (yE<Y(li)),
  48.         V(hi,1:n) = E;       % Replace a vertex:
  49.         Y(hi) = yE;
  50.       else
  51.         V(hi,1:n) = R;       % Replace a vertex:
  52.         Y(hi) = yR;
  53.       end
  54.     end
  55.   else
  56.     if (yR<Y(hi)),
  57.       V(hi,1:n) = R;         % Replace a vertex:
  58.       Y(hi) = yR;
  59.     end
  60.     C = (V(hi,1:n)+M)/2;     % Construct vertex C:
  61.     yC = feval(Fn,C);
  62.     C2 = (M+R)/2;            % Construct vertex C2:
  63.     yC2 = feval(Fn,C2);
  64.     if (yC2<yC),
  65.       C = C2;                % Replace a vertex:
  66.       yC = yC2;
  67.     end
  68.     if (yC<Y(hi)),
  69.       V(hi,1:n) = C;         % Replace a vertex:
  70.       Y(hi) = yC;
  71.     else
  72.       for j = 1:n+1,         % Shrink the simplex:
  73.         if (j~=lo),
  74.           V(j,1:n) = (V(j,1:n)+V(lo,1:n))/2;
  75.           Z = V(j,1:n);
  76.           Y(j) = feval(Fn,Z);
  77.         end
  78.       end                    % End of Shrink.
  79.     end
  80.   end                        % End of Improve.
  81.   [mm lo] = min(Y);          % Order the vertices:
  82.   [mm hi] = max(Y);
  83.   li = hi;
  84.   ho = lo;
  85.   for j = 1:n+1,
  86.     if (j~=lo & j~=hi & Y(j)<=Y(li)), li=j; end
  87.     if (j~=hi & j~=lo & Y(j)>=Y(ho)), ho=j; end
  88.   end                        % End of Order.
  89.   cnt = cnt+1;
  90.   P(cnt,:) = V(lo,:);
  91.   Q(cnt) = Y(lo);
  92.   home; format long;         % Print iteration and plot 2-dim case.
  93.   if cnt==1,
  94.     diary output,...
  95.     disp('The results for a Nelder-Mead search.'),...
  96.     disp('     p                  q                  f(p,q)'),...
  97.     diary off,clc; 
  98.   end; 
  99.   if show==1,
  100.     Mx1 = 'Nelder-Mead search iteration No. ';
  101.     Mx2 = '     p                  q                  f(p,q)';
  102.     disp([Mx1,int2str(cnt)]),disp(Mx2),...
  103.     diary output,disp([V(lo,:),Y(lo)]),diary off;
  104.     XS = V(1:n+1,1)'; XSL = [XS,XS(1)];
  105.     YS = V(1:n+1,2)'; YSL = [YS,YS(1)];
  106.     plot(XS,YS,'or',XSL,YSL,'-g');
  107.   end;
  108. end                          % End of the main while loop.
  109. hold off;
  110. snorm = 0;                   % Determine the size of the simplex:
  111. for j = 1:n+1,
  112.   s = norm(V(j)-V(lo));
  113.   if (s>=snorm), snorm = s; end
  114. end
  115. Q = Q';
  116. V0 = V(lo,1:n);
  117. y0 = Y(lo);
  118. dV = snorm;
  119. dy = abs(Y(hi)-Y(lo));
  120.  
  121.